library(tidyverse)
library(plyr)
library(ggthemes)
library(foreign)

setwd("~/Replication Code/")

##---Loading in data ----
anes <- read.dta("./data/anes_subset.dta")

anes16 <- read.dta("./data/anes16_subset.dta")

# Subsetting and stacking data
dat <- subset(anes, select = c("year", "rr_sc", "pid3", "rep", "pid7_sc"))
dat16 <- subset(anes16, f2f == 1, select = c("rr_sc", "pid3", "rep", "pid7_sc"))
dat16$year <- 2016
cdat <- rbind(dat, dat16)


### Figure 1: Scale Means ----
# Generating means
d <- ddply(cdat, 
           .(year, rep), 
           summarise, 
           mean.rr = mean(rr_sc, na.rm = T))

d <- subset(d, !is.na(rep))
d$Party <- factor(d$rep, labels = c("Democrats", "Republicans"))

## Adding correlation to data frame
YEAR <- unique(d$year)
dcorr <- data.frame(year = YEAR, Corr = rep(NA, length(YEAR)))
for(i in 1:nrow(dcorr)){
  df <- subset(cdat, year == YEAR[i])
  if("TRUE" %in% names(table(!is.na(df$rr_sc)))){
    dcorr$Corr[i] <- cor(cbind(df$rr_sc,df$pid7_sc), use = "complete.obs")[2,1]
  }
}
d <- merge(d, dcorr, by = "year")

# Plotting
ggplot(d, aes(x = year, y = mean.rr, color = Party)) +
  geom_line(data = d[!is.na(d$mean.rr),], aes(lty = Party), size = 1.25) +
  geom_point(data = d[!is.na(d$mean.rr),], aes(shape = Party), 
             size = 3) +
  geom_line(data = d[!is.na(d$Corr),], aes(y = Corr), lty = 4,
            color = "black", size = 1.25) +
  theme_hc() +
  scale_color_manual(values = c("black", "grey")) +
  scale_shape_manual(values = c(16, 15)) +
  annotate("text", x = 2002, y = .75,
           label = "Republicans", color = "grey", size = 6) +
  annotate("text", x = 2003, y = .50,
           label = "Democrats", color = "black", size = 6) +
  annotate("text", x = 1998, y = .28,
           label = "Correlation", size = 6) +
  labs(#title = "Racial Resentment by Party",
    #subtitle = "0-1 Scale (Low-High). American National Election Studies Data.\nFace-to-Face Interviews with Non-Hispanic Whites.",
    y = "Scale Mean",
    x = "Year") +
  theme(legend.position = "none",
        axis.title = element_text(size = 20),
        axis.text = element_text(size = 18)) +
  scale_x_continuous(breaks = seq(1986, 2016, 6)) +
  scale_y_continuous(breaks = seq(0,1,.2), limits = c(0,1),
                     sec.axis = sec_axis(~.*1, name = "Correlation",
                                         breaks = seq(0,1,.2)))
# ggsave("./figures/rr_pid_f2f.pdf", height = 6, width = 9.5)

### Figure 2: Distributions ----
## Generating Counts
counts88 <- prop.table(table(cdat$rep[which(cdat$year == 1988 & cdat$rep == 0
                                            | cdat$year == 1988 & cdat$rep == 1)],
                             cdat$rr_sc[which(cdat$year == 1988 & cdat$rep == 0
                                              | cdat$year == 1988 & cdat$rep == 1)]),1)
counts92 <- prop.table(table(cdat$rep[which(cdat$year == 1992 & cdat$rep == 0
                                            | cdat$year == 1992 & cdat$rep == 1)],
                             cdat$rr_sc[which(cdat$year == 1992 & cdat$rep == 0
                                              | cdat$year == 1992 & cdat$rep == 1)]),1)
counts00 <- prop.table(table(cdat$rep[which(cdat$year == 2000 & cdat$rep == 0
                                            | cdat$year == 2000 & cdat$rep == 1)],
                             cdat$rr_sc[which(cdat$year == 2000 & cdat$rep == 0
                                              | cdat$year == 2000 & cdat$rep == 1)]),1)
counts04 <- prop.table(table(cdat$rep[which(cdat$year == 2004 & cdat$rep == 0
                                            | cdat$year == 2004 & cdat$rep == 1)],
                             cdat$rr_sc[which(cdat$year == 2004 & cdat$rep == 0
                                              | cdat$year == 2004 & cdat$rep == 1)]),1)
counts08 <- prop.table(table(cdat$rep[which(cdat$year == 2008 & cdat$rep == 0
                                            | cdat$year == 2008 & cdat$rep == 1)],
                             cdat$rr_sc[which(cdat$year == 2008 & cdat$rep == 0
                                              | cdat$year == 2008 & cdat$rep == 1)]),1)
counts12 <- prop.table(table(cdat$rep[which(cdat$year == 2012 & cdat$rep == 0
                                            | cdat$year == 2012 & cdat$rep == 1)],
                             cdat$rr_sc[which(cdat$year == 2012 & cdat$rep == 0
                                              | cdat$year == 2012 & cdat$rep == 1)]),1)

counts16 <- prop.table(table(cdat$rep[which(cdat$year == 2016 & cdat$rep == 0
                                            | cdat$year == 2016 & cdat$rep == 1)],
                             cdat$rr_sc[which(cdat$year == 2016 & cdat$rep == 0
                                              | cdat$year == 2016 & cdat$rep == 1)]),1)
## Plotting
COL <- c("black", "white")

# pdf("./Figures/rr_dist16_6panel.pdf", width = 8.25, height = 6.75) # change red to white for paper use
#par(mfrow = c(2,2))
par(mfrow = c(2,3))
barplot(counts88, beside = T, col = COL, 
        axes = F, ylim = c(0,.16),
        main = "1988", xlab = "", ylab = "Proportion")
legend("topleft", c("Democrats", "Republicans"), fill = COL, bty = "n")
axis(2, las = 1)
barplot(counts92, beside = T, col = COL, 
        axes = F, ylim = c(0,.16),
        main = "1992", xlab = "", ylab = "Proportion")
axis(2, las = 1)
barplot(counts04, beside = T, col = COL, 
        axes = F, ylim = c(0,.16),
        main = "2000", xlab = "", ylab = "Proportion")
axis(2, las = 1)
barplot(counts04, beside = T, col = COL, 
        axes = F, ylim = c(0,.16),
        main = "2004", xlab = "", ylab = "Proportion")
axis(2, las = 1)
barplot(counts12, beside = T, col = COL, 
        axes = F, ylim = c(0,.16),
        main = "2012", xlab = "", ylab = "Proportion")
axis(2, las = 1)
barplot(counts16, beside = T, col = COL, 
        axes = F, ylim = c(0,.16),
        main = "2016", xlab = "", ylab = "Proportion")
axis(2, las = 1)
# arrows(x0 = , y0 = -0.02, xpd = T, lwd = 2)
# dev.off()


## Testing difference in distributions
Y <- 1988
summary(table(cdat$rep[which(cdat$year == Y & cdat$rep == 0
                             | cdat$year == Y & cdat$rep == 1)],
              cdat$rr_sc[which(cdat$year == Y & cdat$rep == 0
                               | cdat$year == Y & cdat$rep == 1)]))
